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ABSTRACT 



Aims. We establish a criterion for the stability of planetary orbits in stellar binary systems by using Lyapunov exponents 
and power spectra for the special case of the circular restricted 3-body problem (CR3BP). The criterion augments our 
py^ earlier results given in the two previous papers of this series where stability criteria have been developed based on the 

j-f^ Jacobi constant and the hodograph method. 

Methods. The centerpiece of our method is the concept of Lyapunov exponents, which are incorporated into the analysis 
of orbital stability by integrating the Jacobian of the CR3BP and orthogonalizing the tangent vectors via a well- 
ed i' established algorithm originally developed by Wolf et al. The criterion for orbital stability based on the Lyapunov 
exponents is independently verified by using power spectra. The obtained results are compared to results presented in 
the two previous papers of this series. 

Results. It is shown that the maximum Lyapunov exponent can be used as an indicator for chaotic behaviour of planetary 
orbits, which is consistent with previous applications of this method, particularly studies for the Solar System. The 
chaotic behaviour corresponds to either orbital stability or instability, and it depends solely on the mass ratio /i of the 
binary components and the initial distance ratio po of the planet relative to the stellar separation distance. Detailed 
case studies are presented for /i = 0.3 and 0.5. The stability limits are characterized based on the value of the maximum 
^ • Lyapunov exponent. However, chaos theory as well as the concept of Lyapunov time prevents us from predicting exactly 

' when the planet is ejected. Our method is also able to indicate evidence of quasi-periodicity. 

Conclusions. For different mass ratios of the stellar components, we are able to characterize stability limits for the 
CR3BP based on the value of the maximum Lyapunov exponent. This theoretical result allows us to link the study 
■ of planetary orbital stability to chaos theory noting that there is a large array of literature on the properties and 

significance of Lyapunov exponents. Although our results are given for the special case of the CR3BP, we expect that 
it may be possible to augment the proposed Lyapunov exponent criterion to studies of planets in generalized stellar 
binary systems, which is strongly motivated by existing observational results as well as results expected from ongoing 
and future planet search missions. 
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5— i ' 1. Introduction The study of planetary orbital stability is timely for 

various astronomical reasons. First, although most plan- 
A classical problem within the realm of orbital stability ets are found in wide binaries, various cases of plan- 
studies for theoretical and observed planets in stellar binary ets in binaries with separation dis tances of less than 
systems is the circular restricted 3-body problem (CR3BP) 30 AU have also been id e ntified (e.g.. iPatience et al.l 20021: 
(e.g., ISzebehelvl Il967t iRovl 120051) . The CR3BP describes lEggenberger et al.l 12004 lEggenberger fc Udrv 1 12010, and 



the motion of a body of negligible^ mass moving in the references therein). Second, many more cases of planets in 

gravitational field of the two massive primaries considered stellar binary systems are expected to be discovered not- 

here to be two stars. The stars move in circular orbits about ing that binary (and possibly multiple) stellar systems oc- 

the center of mass and their motion is not influenced by the cur in high frequency i n the local Galact i c neighbourhood 

third body, the planet. Furthermore, the initial velocity of (iDuauennov fc Mavorl Il99lt lLadal 120061: iRaghavan et al.l 

the planet is assumed in the same direction as the orbital |2006l ). Moreover, orbital stability studies of planets around 

velocity of its host star, which is usually the more massive stars, including binary systems, are highly significant in 

of the two stars. consi deration of ongoing and future planet search missions 



consider ation of on going and future planet search 
(e.g., ICatanzarite et alj|2006t ICockell et a"L1l2009f) . 



Send offprint requests to: Z. E. Musielak The CR3 B P has been s tudied in de tail by 

Negligible mass means that although the body's motion is |Dvorald (11984T). iDvorakl (119861) iRabl fc Dvorak! (Il988t). 



influenced by the gravity of the two massive primaries, its mass ISmith fc Szebehefvl (Il993h ~ iPilat-Lohinger fc Dvorak! 
is too low to affect the motions of the primaries. (2002|), and [Mardling (2007]). It has also been the focus of 
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the pr evious papers in this series. In Paper I, lEberle et al.l 
(|2008f) obtained the planetary stability limits using a 
criterion based on Jacobi's integral and Jacobi's constant. 
The metho d, also related to the concept of Hill stability 
(|Rovll2005t ). showed that orbital stability can be guaran- 
teed only if the initial position of the planet lies within 
a well-defined limit determined by the mass ratio of the 
stellar components. Additionally, the stability criterion 
was found to be also related to the borders of the "zero 
velocity contour" (ZVC) and its topology assessed by using 
a synodic coordinate system regarding the two stellar 
components. 

In Paper II, lEberle fc Cvmtel (|2010t ) followed another 
theoretical concept based on the hodograph eccentricity cri- 
terion. This method relies on an approach given by differ- 
ential geometry that analyzes the curvature of the plan- 
etary orbit in the synodic coordinate system. The cen- 
terpiece of this method is the evaluation of the effective 
time-dependent eccentricity of the orbit based on the hodo- 
graph in rotating coordinates as well as the calculation of 
the mean and median values of the eccentricity distribu- 
tion. This approach has been successful in mapping quasi- 
periodicity and multi-periodicity for planets in binary sys- 
tems. It has also been test ed by comparing its theor etical 
predictions with wor k by iHolman fc Wiegertl ([1999) and 
iMusielak et al.l (|2005T) in regard to the extent of the region 
of orbital stability in binary systems of different mass ra- 
tios. 

Previously, the work by lEberle fc Cunt3 <|2010h dealt 
with the assessment of short-time orbital stability, encom- 
passing time scales of 10 3 yrs or less. One of the findings 
was the identification of a quasi-periodic region in the stel- 
lar mass and distance (p, po) parameter space (see Sect. 
2.1). Due to the relatively short time scales considered in 
the previous study, there is a strong motivation to revisit 
the onset of orbital instability using longer time scales and 
different types of methods. The focus of this paper is the 
analysis of orbital stability by Lyapunov exponents, which 
are among the most commonly used numerical tools for 
investigati ng chaotic beh aviour of different dynamical sys- 
tems (e.g.. iHilbornl 11993 ). The exponents have already re- 
peatedly bee n used in orbita l mechanics studies of the Solar 
System ( e.g.. lLissauerlll999t iMurrav fc Hormanll200lD . For 
example, iLissauerl (|1999[ ) discussed the long-term stability 
of the eight Solar System planets, while also considering 
previous studies. He concluded that the Solar System is 
most likely astronomically stable, in the view of the limited 
life time of the Sun; however, the orbits of Pluto and of 
many asteroids may become unstable soon after the Sun 
becomes a white dwarf. 

The numerical procedure of c omputing the Ly apunov 
exponents has been developed by IWolf et al.l (|1985( ) based 
on previous work bv lBenettin et all ( 1980t ). Hence, the main 
objective of the present paper is to establish the Lyapunov 
exponent criterion an d verify it by performing the power 
spectra analysis (e.g.. lHilborn| [T994) as well as by compar- 
ing the obtained results to those presented in Paper I and II. 
Our newly established criterion is then used to investigate 
the stability of planetary orbits in stellar binary systems 
(approximated here as the CR3BP) of different mass and 
distance ratios. The methodology of our study of orbital 
stability based on the Lyapunov exponent criterion is pre- 
sented and discussed in §2. Detailed model simulations are 
described in §3, and our conclusions are given in §4. 



2. Theoretical approach 

2.1. Basic equations 

In th e following, we consider th e so-called coplanar CR3BP 
fe.g.. lSzebehelvill967l:lRovll2005l) . which we will define as fol- 
lows. Two stars are in circular motion about their common 
center of mass and their masses are much larger than the 
mass of the planet. In our case, the planetary mass is as- 
sumed to be 1 x 1CP 6 of the mass of the star it orbits; also 
note that the planetary motion is constrained to the orbital 
plane of the two stars. Moreover, it is assumed that the ini- 
tial velocity of the planet is in the same direction as the 
orbital velocity of its host star, which is the more massive 
of the two stars, and that the starting position of the planet 
is to the right of the primary star (3 o'clock position), along 
the line joining the binary components. 

Fo llowing the conventions described by lEberle et al.l 
(2008), we write the equations of motion in terms of the 
parameters p and po with p and a = 1 — p being related to 
the ratio of the two stellar masses m\ and mi (see below). 
Moreover, i?o denotes the initial distance of the planet from 
its host star, the more massive of the two stars with mass 
mi, whereas D denotes the distance between the two stars, 
allowing us to define pq. In addition, we use a rotating refer- 
ence frame, which also gives rise to Coriolis and centrifugal 
forces. The equations of motion are given as 



x 

y 

z 



u 
v 
w 

2v 



= -2u 



— a — „ p— 

y y 

' 1 '2 

Z 



where 



TO2 



a 
r\ 
r\ 

Po 



mi + 7712 

l-p 

(x - p) 2 + y 2 + z 2 

(x + a) 2 + y 2 + z 2 

Ro 
D 



(1) 
(2) 
(3) 

(4) 
(5) 
(6) 

(7) 

(8) 
(9) 
(10) 

(11) 



The variables in the above equations describe the po- 
sition of the planet, which in essence constitutes a test 
particle. Its position is defined in Cartesian coordinates 
{x, y, z}. We denote the time derivative or velocity of a 
coordinate using the dot notation {i = 4|} . We also rep- 
resent the set of second order differential equations, the 
equations of motion, by a set of first order differential equa- 
tions. The velocity is defined by the coordinates {w, i>, w} 
whose time derivatives are the accelerations. By defining 
the mass ratio and using normalized coordinates, we can 
define the distances {ri, ri} with reference to the location 
of the stars in the rotating coordinate system. 

We enumerate the co-linear Lagrange points in the syn- 
odic frame by the order of which the ZVC opens. The point 
between the stars opens first; therefore, we denote it as LI. 
The point to the left of the star that does not host the 



Quarles et al.: The instability transition for the restricted 3-body problem. III. 



3 



planet opens second; thus, we denote it L2. The point to 
the right of the star hosting the planet opens third and it is 
denoted L3. The two Trojan Lagrange points which are of 
lesser importance to our study are L4 and L5 (see Fig. [IJ. 

2.2. Lyapunov exponents 

A fundamental difference between stable and unstable plan- 
etary orbits is that two nearby trajectories in phase space 
will diverge as a power law (usually linear) for the former 
and exponentially for the latter. The parameter that is used 
to measure this rate of divergence is ca lled Lyapunov ex - 
ponent as it was originally intro duced bv lLvapunovl (|l907f k 
see also lBaker fc Gollub I (|1990f ). A dynamical system with 
n degrees of freedom is represented in 2n phase space; 
thus, to fully determine the stability of the system all 2n 
Lyapunov exponents must be calculated. The Lyapunov 
exponents are the most commonly used tools to deter- 
mine the onset of chaos and chaotic beh aviour of both 
dissipative (e.g.. lMusielak fc Musielakll2009l and references 
therein) and non-dissip ative systems of orbital mechan- 
ics (e.g.. lLissauerlll99"9T and references therein). The posi- 
tive Lyapunov exponents measure the rate of divergence of 
neighbouring orbits, whereas negative exponents measure 
the convergence rates between stable manifolds. For dissi- 
pative dynamical syste ms the sum of all Lyapunov expo- 
nents is less than (e.g.. lMusielak fc Musielakll2009() ; how- 
ever, for Hamilt onian (non-dis sipative) systems the sum is 
equal to (e.g.. lHilbornlll994ft . 

Specific applications of the Lyapunov exponents to 
the circular restricted 3-bod y problem were d is cusse d 
by many authors, i ncluding [G onczi & Froesch ll (119811). 
iJefferys fc Yl (ll983ri. ITecar et al.l (fl992h. iMilani fc Nobili 
dl992 | ) . ISmith fc Szebehelvi (|l993h . iMurrav fc Holman 
(|2001[) . and others. Some of these authors also considered 
the so-called Lyapunov time, which measures the e-folding 
time for the divergence of nearby trajectories. It should 
be noted that the Lyapunov time is not well-defined for 
cases when close encounters between a planet and one of 
the s tars occur or w hen the planet is ejected from the sys- 
tem (|Lissauerlll999f) . We shall treat such cases with special 
caution in this paper. 

The previously obtained results clearly show that the 
Lyapunov exponents can be calculated for the case of the 
CR3BP considered in this paper, for which we have 2n = 6. 
This requires a state vector for the system containing 6 el- 
ements. Details about the nature of the state vector are 
given in Appendix A. From the equations of motion, the 
Jacobian J can be determined, which will be the founda- 
tion of how the Lyapunov exponents will be determined. In 
addition, the nature of this Jacobian will elude to certain 
properties of the Lyapunov exponents. 

For Hamiltonian systems the trace of the Jacobian 
should equal zero, Tr J = 0. This requires that the CR3BP 
should have either all zero diagonal elements or an even 
amount of 3 positive/negative diagonal elements. As a re- 
sult this forces the sum of the Lyapunov exponents to be 
also zero. Then we should expect some symmetry in the 
spectrum of the Lyapunov exponents. Since this is indeed 
the case, the Lyapunov exponent spectra shown in this 
paper will only present the positive Lyapunov exponents 
and omit the negative Lyapunov exponents as they do not 
give any additional information about the system. It is the 
convention that positive Lyapunov exponents indicate that 
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Fig. 1. Locations of the Lagrange points LI, L2, L3, L4, 
and L5 as used in the CR3BP. The label P denotes the 
primary star and the label S denotes the secondary star. 



both dissipative (e.g., Hilbornl 1 1994T) and non-dissipative 



(e.g., lOzorio de Almeida 19901 ) systems are chaotic. Here 



the end value of the Lyapunov exponents will be used to 
make the distinction between chaotic and non-chaotic or- 
bits. 

It is well known that the largest Lyapunov exponent 
is sufficient t o determine the m agnitude of the chaos in a 
system (e.g., IWolf et al.|[l985l) . Therefore, this study will 
consider only the effect of the maximum Lyapunov expo- 
nent because it will show the greatest effect of chaos on the 
system. This motivates us to also examine how large the 
maximum Lyapunov exponent can be before introducing 
enough chaos to affect the stability of the CR3BP. 

3. Results and discussion 

We performed simulations for stellar mass ratios from p = 
0.0 to 0.5 in increments of 0.01. A Yoshida s ixth order sym - 
plectic integration scheme was used (e.g.. lYoshidalll990D . 
As a measure of the precision of the integration scheme, 
we note that a time step of e = 10~ 4 yrs is used for the 
individual cases. However, for producing plots of the entire 
parameter space we adopt a time step of e = 10~ 3 yrs as 
the smaller time step did not noticeably enhance the qual- 
ity of the plots. The order of error in energy for each time 
step was 10 -14 and 10 -10 , respectively. We performed sim- 
ulations for different time limits, which range from 10 to 
10 5 yrs in increments of powers of 10. We also performed 
case studies with time limits of 10 6 and 10 7 yrs, although we 
do not expect significant changes to occur at those longer 
time scales. By using different time scales we can estimate 
when certain phenomena occur and ascertain how they will 
affect other regions over longer periods of time. 

We display runs for selected initial conditions (i.e., start- 
ing distances po of the planet from the stellar primary) 
corresponding to p = 0.3 in Figs. [5] to 0] and to p — 0.5 in 
Figs. [5] to [7] In Figs. [5] to El we first present the planetary 
orbits in the synodic coordinate system (X*, Y*). Secondly 
in each figure, the Lyapunov spectra of the simulations are 
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Table 1. Errors in the Jacobi constant (JCE) for models of p — 0.3 at different time intervals. 



Po 


JCE(J 


) 


JCE(5 


) 


JCE(2f) 


JCE(r 


) 


0.355 
0.474 
0.595 


1.0740E- 
1.4498E- 
8.9819E- 


10 
10 
13 


2.4747E- 
9.7917E- 
5.2923E- 


10 
11 
13 


1.2403E-10 
3.7723E-11 
5.2001E-13 


1.8824E- 
1.0067E- 
7.1001E- 


10 
11 
06 



Note: t denotes the runtime of the simulation. We show the stable cases of po = 0.355 and 0.474 with r = 10 5 binary orbits as 
well as the unstable case of po = 0.595 with r = 150.64 binary orbits. 



Table 2. Maximum Lyapunov exponent study for the models of fi = 0.3 at different time intervals. 



PO 


MLE (10 


2 ) 


MLE (10 


3 ) 


MLE (10 


4 ) 


Cmcdian 


ZVC 


0.20 


1.1099E- 


•1 


1.1493E- 


2 


1.1446E- 


3 


0.027 




0.30 


1.0322E- 


1 


1.0103E- 


2 


1.0691E- 


3 


0.17 




0.40 


9.8084E- 


2 


1.0050E- 


2 


6.8986E- 


4 


0.58 


LI 


0.474 


7.7766E- 


2 


9.4232E- 


3 


9.2624E- 


4 


0.75 


LI, L2 


0.50 


8.5390E- 


2 


9.2048E- 


3 


7.9830E- 


4 


0.85 


LI, L2 


0.595 


1.5793E- 


1 


1.4175E- 


1* 






1.14 


LI, L2, L3 


0.60 


8.9533E- 


1* 










1.10 


LI, L2, L3 



Note: Elements without data indicate simulations that ended due to the energy criterion, thus representing planetary catas- 
trophes. Elements with an asterisk (*) indicate that the simulation ended before the allotted time. Also, e me di an represents the 
median of the eccentricity distribution obtained for the curvature of the planetary orbits in the synodic coordinate system (see 
Paper II). The last column indicates the Lagrange point(s) where the zero-velocity contour (ZVC) is open (see Paper I). 



shown using a logarithmic scale for the Y*-axis. Lastly in 
each figure, the time series power spectra of the simulations 
are shown with their normalized amplitudes. The power 
spectra were obtained through the usage of a FFT subrou- 
tine in Matlab® that uses the X*-component separation 
distance as a function of time and converts the output fre- 
quencies to periods. The selected starting distance ratios 
po of the planet are 0.355, 0.474, 0.595 for p = 0.3 and 
0.290, 0.370, 0.400 for /i = 0.5. Note that p indicates the 
relative initial distance of the planet given as Rq/D, with 
D = 1 AU as the distance between the two stars and Rq 
as the initial distance of the planet from its host star, the 
primary star of the binary system. 

Using the method of Lyapunov exponents , we are able 
to ver ify a nd extend the methods d escribed bv lEberle et al.l 
(|2008l) and lEberle fc Cuntzl (|2010D . Absolute orbital stabil- 
ity can be more rigorously shown through the Lyapunov 

exponent method. It occurs when po < p^\ where p^ rep- 
resents the point at which the initial condition results in a 
ZVC that opens at LI. For larger values of p n, stability is 
not gu aranteed due the behaviour discussed bv lEberle et al.l 
(2008(); this result is also consistent with our Lyapunov ex- 
ponent criterion for orbital stability. 

In this paper, the main criterion for orbital stability is 
the Lyapunov exponent criterion, which is based on the 
maximum Lyapunov exponent. From a theoretical point of 
view, an orbit is stable when all Lyapunov exponents are 
exactly zero. Obviously, this 'perfect' criterion for orbital 
stability will be very difficult to reach numerically because 
it would require an infinite simulation time. Based on our 
finite simulation times, we obtain 3 positive and 3 negative 
Lyapunov exponents, and their sum becomes close to zero 
within the limits of our numerical simulations. Hence, in or- 



der to determine orbital stability numerically, we must look 
at the values of the three positive Lyapunov exponents at 
the beginning and at the end of our simulations. By com- 
paring these values, we determine whether the exponents 
decrease in time, and if so what is the rate of their decrease, 
or whether they stay approximately constant in time. 

Using this information, we are able to identify the maxi- 
mum Lyapunov exponent and adopt it as our primary indi- 
cator of orbital stability. We classify an orbit as stable if the 
initial value of the maximum Lyapunov exponent is below 
a certain threshold and if it decreases at a 'reasonable rate' 
(see below); otherwise, the orbit is classified as unstable. 
In our plots of Lyapunov spectra (see the second panels in 
Figs. [5] to [7]), we present all three positive Lyapunov expo- 
nents to show their values and study how they change in 
time. 

In addition to the Lyapunov exponent criterion for or- 
bital stability, we also use the so-called orbital energy cri- 
terion, which requires that the kinetic energy should not 
exceed twice the value of the potential energy. This is eval- 
uated in our numerical simulations during each time step. A 
failure of this criterion would imply a break in c o nserv ation 
of the Jacobi constant as detailed by ISzebehelvl (|1967|) and 
shown numerically in Table [TJ It should also be noted that 
the cases presented in Table Q] are also presented in Figs. [5] 
to [4] There are two important points associated with this 
criterion. First, our systems are Hamiltonian, which means 
that their energy must be conserved. On the other hand, 
the fact that the energy conservation of the planet may be 
broken because we are neglecting the effects of the third 
mass on t h e larg er masses has already been discussed by 
ISzebehelvl (jl967l ). Second, even if the energy is only ap- 
proximately conserved in our numerical simulations, we still 
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Table 3. Maximum Lyapunov exponent study for the models of /i = 0.5 at different time intervals. 



Po 


MLE (10 2 ) 


MLE (10 3 ) 


MLE (10 4 ) 


^median 


ZVC 


0.25 


8.7189E-2 


1.1353E-2 


1.0009E-3 


0.20 




0.29 


1.0647E-1 


9.5469E-3 


9.9141E-4 


0.45 


LI 


0.30 


7.6275E+0* 






0.67 


LI 


0.35 


5.2937E+0* 






0.82 


LI 


0.37 


6.1946E+0* 






0.98 


LI 


0.40 


1.0675E-1 


9.6613E-3 


8.8900E-4 


0.71 


LI 


0.43 


2.8614E-1 


2.9420E-1* 




0.87 


LI 


0.50 


7.9684E-1* 






0.089 


LI, L2, L3 



Note: Elements without data indicate simulations that ended due to the energy criterion, thus representing planetary catas- 
trophes. Elements with an asterisk (*) indicate that the simulation ended before the allotted time. Also, e me di an represents the 
median of the eccentricity distribution obtained for the curvature of the planetary orbits in the synodic coordinate system (see 
Paper II). The last column indicates the Lagrange point (s) where the zero- velocity contour (ZVC) is open (see Paper I). 



request that the Jacobi constant remains constant, which is 
a required stability condition for CR3BP. To be consistent 
with this criterion (see also Paper I), we stop our numerical 
simulation once changes (even small) in the Jacobi constant 
occur; such cases are depicted in Tables [2] and [3l 

Alike in the previous paper bv lEberle et al.l ()2008[ ). we 
are able to identify the same regions of orbital stability, in- 
stability, and domains of quasi-periodic orbital stability. A 
key difference, however, is that the orbit diagrams depicted 
in Figs. [2] to [7] have been simulated for 10 5 years, whereas 
the corresponding figures in the previous paper have been 
simulated for only 10 3 years. The power spectra are shown 
to indicate how we determined the region of quasi-periodic 
orbital stability including the correct magnitude. 

We begin by classifying Figs. [2 [3j [5j and [Jj as stable 
candidate configurations. Three of the cases show similar 
behaviour in the maximum Lyapunov exponent as shown 
in Figs.[3j},[5)D, and[7j3. These Lyapunov spectra have a com- 
mon trend by starting at a maximum value on the order of 
10 _1 and converging, albeit slowly, to smaller orders of ten. 
Figure [3] shows a somewhat different behaviour compared 
to the other cases in its class of stable candidates. This 
case establishes a quasi-periodic orbit, which illustrates a 
(3:1) resonance in the power spectrum and reveals the 
same trend of stability in the Lyapunov exponent spectrum. 
Figure shows an additional variance in behaviour due to 
the elevated nature of the maximum Lyapunov exponent. 
This shows that a limited amount of chaos exists in the 
system while remaining stable for the full simulation time. 
This case also demonstrates a (4 : 1) resonance in the power 
spectrum. 

In contrast, Fig. [4] demonstrates a case of instability. 
The Lyapunov spectrum shows a different nature than that 
of the preceding cases. In Fig. 0t> it begins at a maximum 
value greater than 1 and converges to a value between 10 _1 
and 1. By establishing the preceding cases as stable cases, 
we can contrast the final values of the corresponding maxi- 
mum Lyapunov exponents. In the unstable case of Fig. \4)p, 
it is two orders of magnitude greater than in the stable 
cases of Figs. [5b, [5b, and [7b. However, in comparison with 
Fig-Hb we nn d only a difference of one order of magnitude. 

Figure [6] has been examined in a similar manner. Having 
already classified Fig. [5] and [7j as stable configurations, we 
emphasize that they reveal similar trends in the maximum 



Lyapunov exponent as given by Fig. [3] However, Fig. [5] 
shows a different orbit diagram and is described by a max- 
imum Lyapunov exponent that gives the outcome of insta- 
bility. This case conveys a much noisier power spectrum 
along with a maximum Lyapunov exponent on the order 
of 1 or greater. Finer detail is illustrated in Table [TJ and 
Tabled Particularly Table[5]shows a boundary in the maxi- 
mum Lyapunov exponent where values near 0.1 and smaller 
for the first 100 years tend toward stability. 

Figure [5] conveys the bigger picture for the overall pa- 
rameter space. It represents contour plots of the maximum 
Lyapunov exponent with respect to the (p, po) parame- 
ter space in linear and logarithmic scale. The crosses de- 
pict initial conditions for runs that prematurely ended due 
to the orbital energy criterion. This demonstrates where 
the regions of instability occur as reflected by the respec- 
tive colour coded scale of the maximum Lyapunov expo- 
nent. Comparing F ig. [SJ (left) to the previous result by 
lEberle et all ([20081 ) . it can be seen that the main regions 
of stability remain the same. Other noteworthy features in- 
clude that the instability islands present near po — 0.4 have 
grown as the simulation time has been increased as well as 
the existence of a plateau near po = 0.48. The colour scale 
in Fig. [5] (left) has a maximum colour of dark red at a max- 
imum Lyapunov exponent of 0.15. Some regions appear to 
be coloured black; however, this does not correspond to the 
adopted colour scale as it is caused by the finite contour line 
width due to the close proximity of the lines. Therefore, the 
plateau has a peak value of a maximum Lyapunov exponent 
near 0.15. Considering the diagrams in Fig. [31 we conclude 
that a region of quasi-periodicity exists on this plateau. 

Figure [S] (right), the contour plot in logarithmic scale, 
is shown to display some of the finer details pertaining 
to structure inside the stability regions. The blue-green 
coloured regions demonstrate areas of possible stable chaos 
that are hidden in Fig. [5] (left). Some other smaller contours 
are also produced in the stability region, but the colour cod- 
ing in general shows little change, which is partially due to 
the chosen spacing between contours. The average differ- 
ence in height between these levels is —0.15; on the loga- 
rithmic scale, it indicates a change by a factor of 1.41 in 
the maximum Lyapunov exponent. 

Inspecting Fig. [8] also allows a comparison with previous 
results, particularly results of Paper I. There it was shown 
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that absolute orbital stability occurs when po < p a , where 

Pq 1 '' represents the point at which due to the initial condi- 
tion the ZVC opens at LI. This condition is represented by 
the line of Ci in Fig. [5] It is evident that the Lyapunov 
exponent criterion is almost perfectly consistent with this 
condition. Additionally, as discussed in Paper I and ref- 

■ (3) 

erences therein, no stability can be expected if po > p Q 
because now the ZVC is open at LI, L2, and L3. This find- 
ing is also almost perfectly reflected through the behaviour 
of the maximum Lyapunov exponent as a orbital stability 
criterion. Both panels of Fig.[8]show no stability regions on 
the right side of C3, which is reflective of the opening of the 
ZVC at Lagrange points LI, L2, and L3. 



4. Conclusions 

We present a detailed case study of the CR3BP with anal- 
yses through the usage of the Lyapunov exponents. We are 
able to characterize stability limits based on the value of 
the maximum Lyapunov exponent at the end of each simu- 
lation. Cases where the maximum Lyapunov exponent ex- 
ceeds a value of 0.15 indicate that the planet will experience 
an event that causes the orbital velocity to decrease. These 
events include intersecting a Lagrange point or the ZVC. 
After such an event, the planet will experience a series of 
near misses (or collisions) with one of the stars in the binary 
system leading to overall instability. Chaos theory and the 
concept of Lyapunov time prevent us from predicting ex- 
actly when the planet will be ejected. Using the Lyapunov 
time as a measure of the length of predictive time, we can 
show that this relationship is proportional to the inverse of 
the maximum Lyapunov exponent. Using our critical value 
of 0.15, the unstable systems will lose predictability within 
(0.15) -1 = 6| years or less. After that time, it is unknown 
when the planet will be ejected as this could take many 
multiples of the Lyapunov time to occur. 

The method also shows evidence of a region of quasi- 
periodicity. This is a region where the maximum Lyapunov 
exponent remains near the critical value without exceeding 
that value. This is shown for cases simulated for 10 5 years. 
Based on this time scale, we can conclude that our quasi- 
periodic plateau represents a region of stable or quasi-stable 
chaos. This is due to the Lyapunov times being much less 
than the simulation run time. The chaos is shown through 
the value of the maximum Lyapunov exponent, whereas 
the stability is shown through the motion of the planet 
over longer time scales. The lack of near miss events and 
encounters with regions of decreasing velocity prohibits the 
planet from becoming orbitally unstable. 

Comparisons to previously established criteria for sta- 
bility show that our results are consistent with previously 



obtained stability limits (e.g.. iHolman fc Wiegert 
Musielak et"aHl2005HCuntz et alj|2007t lEberle et al.l 



1999: 



2008; 



Eberle fc Cuntl2010l) . In lEberle etail (|2008l ). Paper I, the 
onset of in stability was related to the topology of the ZVC, 
whereas in lEberle fc Cuntzl (|2010f ). Paper II, it was shown 
that the onset of orbital instability occurs when the median 
of the effective eccentricity distribution exceeds unity. Both 
results are consistent with the findings of Paper III, i.e., the 
inspection of the maximum Lyapunov exponent. However, 
the latter offers the advantage to link the study of orbital 
instability for the CR3BP (and any subsequent generaliza- 



tion, if available) to chaos theory, including the evaluation 
of different types of chaos. 

Although our results have been obtained for the special 
case of the CR3BP, we expect that it may also be possi- 
ble to augment our findings to planets in generalized stel- 
lar binary systems. Desired generalizations should include 
studie s of the elliptical restricted 3-bo d y problem (ER3BP) 
(e.g iPilat-Lohinger fc Dvorak] 120021: ISzenkovits fc Mako 
2008) as well as of planets on inclined orbits ([Dvorak et al. 
2007), noting that especially the former have important 
applications to real existing systems in consideration of 
the identified stellar and planetary orbital parameters 
(|Eggenberger fc Udrv II20T0I) . 
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Appendix A 

Basic concepts and definitions 

From stability analysis we can use a Jacobi an matrix to 
show how a state vector x evolves in time (see lTsonislll992l 
for details). The governing equations are 

x = J • x 



\ 9fi 3/„ I 

This differential equation has a general solution, which is 

x (t) — e^ 1 x(0). Assuming J has distinct eigenvalues, we 
can find a matrix U to diagonalize J to form a diagonal 
matrix D such as 

IT 1 J U= D 

By rewriting the above equation we obtain the more useful 
form given as 

J = ud ir 1 

Using the multiplication theorem, we can develop relations 
to obtain the eigenvalues of J yielding 

dct J = (det U) (dct D) (dct IT 1 ) = det D 

= A1A2 • • • A„ 

Furthermore, we find 

Tr J = Tr D = A a + A 2 + . . . + A„ 

This can be generalized from J to a function / (J) such as 
det /(J) = /(A 1 )/(A 2 )-.-/(A„) 
Tr/(J) = /(A 1 ) + /(A 2 ) + --- + /(A„) 

From our general solution, we can assume /(J) = e^ f . 
Therefore, we find 



det e Jt = e (V+A 3 +...+A„)t 



,(TrJ)t 



A volume of perturbations in phase space will be conserved 



if 



det 



,(TrJ)t 



= 1 or Tr J = for Hamiltonian systems. A 



dissipative system will have 



dct e 



(Tr J)( 



< 1 or Tr J < 0. 



The eigenvalues Ai, A2, ■ • • , A„ are the Lyapunov exponents 
of the flow or the logarithms of eigenvalues of J. 
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From these equations we can define 6 dimensional tangent 
vectors x, and their derivatives Xj where i = 1, . . . , 6 as 

Xi {Xi : Ui: %ii Uij Vi: 

Xj = {ui,Vi,Wi,Ui,Vi,Wi} T 

Also we can define a Jacobian matrix, J, given as 



/ 1 0\ 

1 

1 

w w aa o 2 (J 

ox oy oz 

§ % If -2 

\f f I o 00/ 



Wolf method 



Now we describe the scheme by which W olf et ahl £1985) 
determine the Lyapunov spectrum and its application to 
the CR3BP. The Wolf method follows a basic algorithm. 
The first step is to initialize a state vector of 6 ele- 
ments. Then, the tangent vectors need to be initialized to 
some value. We choose to have all tangent vectors to be 
unit vectors for simplicity. This means that the elements 
{xi, 2/2, ^3, U4, t>5, wg} = 1 and all other elements will be 
equal to zero. 

The next step consists of a loop that will use a standard 
integrator akin to the Runge-Kutta schemes to determine 
how the state and tangent vectors will change within a time 
step. This will continue for an adequate number of steps so 
that the tangent vectors become oriented along the flow. 
When this has been accomplished, it is necessary to per- 
form a Gram-Schmidt Renormalization (GSR) to orthogo- 
nalize the tangent space. Thereafter, we take the logarithm 
of the length of each tangent vector to obtain the Lyapunov 
exponents. We then continue the loop. 

Procedure 

We use the first tangent vector, Xi, to define the basis for 
the GSR process. Thus, the first step consists in normalizing 
this vector. With the vectors of the new orthonormal set of 
tangent vectors, denoted by primes (/), we find 



Xn = 



„Xi|| 

_ Xa-Vxa^xAXj 



x 6 -(x 6 , X5)x5-...-(x 6 , Xi)x' t 

|ix 6 -(x 6 , x;)x;-...-(x 6 . x', )x; i 



From this new set of tangent vectors, the Lyapunov ex- 
ponents will be calculated considering the lengths of each 
vector. Therefore, we find 



T 



log- 
where A = x D = 0. 



E 



x * \ Xi ' x J-i/ x j-i 



Baker, G. L., & Gollub, J. P. 1990, Chaotic Dynamics: An 
Introduction (Cambridge: Cambridge University Press) 

Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J.-M. 1980, 
Meccanica, 15, 9 

Catanzarite, J., Shao, M., Tanner, A., Unwin, S., & Yu, J. 2006, PASP, 
118, 1319 

Cockell, C. S., Leger, A., Fridlund, M., et al. 2009, Astrobiology, 9, 1 
Cuntz, M., Eberle, J., & Musielak, Z. E. 2007, ApJ, 669, L105 
Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 
Dvorak, R. 1984, Celest. Mech. Dyn. Astron., 34, 369 
Dvorak, R. 1986, A&A, 167, 379 

Dvorak, R., Schwarz, R., Suli, A., & Kotoulas, T. 2007, MNRAS, 382, 
1324 

Eberle, J., & Cuntz, M. 2010, A&A, 514, A19 (Paper II) 
Eberle, J., Cuntz, M., & Musielak, Z. E. 2008, A&A, 489, 1329 
(Paper I) 

Eggenberger, A., & Udry, S. 2010, in Planets in Binary Star Systems, 

ed. N. Haghighipour (New York: Springer), 19 
Eggenberger, A., Udry, S., & Mayor, M. 2004, A&A, 417, 353 
Fouchard, M., Lega, E., Froeschle, Ch., & Froeschle, CI. 2002, Cel. 

Mech. Dyn. Astron., 83, 205 
Froeschle, C, Lega, E., & Gonczi, R. 1997, Cel. Mech. Dyn. Astron., 

67, 41 

Gonczi, R., & Froeschle, C. 1981, Cel. Mech. Dyn. Astron., 25, 271 
Hilborn, R. C. 1994, Chaos and Nonlinear Dynamics (Oxford: Oxford 

University Press) 
Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 
Jefferys, W. H., & Yi, Z.-H. 1983, Cel. Mech. Dyn. Astron., 30, 85 
Lada, C. J. 2006, ApJ, 640, L63 

Lecar, M., Franklin, F., & Murison, M. 1992, AJ, 104, 1230 
Lissauer, J. J. 1999, Rev. Mod. Phys., 71, 835 

Lyapunov, M. A. 1907, Ann. Fac. Sci., University of Toulouse, 9, 203 
Malhotra, R. 1998, in Solar System Formation and Evolution, ed. D. 

Lazzaro, R. Vieira Martins, S. Ferraz-Mello, & J. Fernandez (San 

Francisco: ASP), 149, 37 
Mardling, R. A. 2007, in Dynamical Evolution of Dense Stellar 

Systems, ed. E. Vesperini, M. Giersz, & A. Sills, Proc. IAU Symp. 

246, 199 

Milani, A., & Nobili, A. M. 1992, Nature, 357, 569 
Murray, N., & Holman, M. 2001, Nature, 410, 773 
Musielak, Z. E., Cuntz, M., Marshall, E. A., & Stuit, T. D. 2005, 

A&A, 434, 355; Erratum, 480, 573 [2008] 
Musielak, Z. E., & Musielak, D. E. 2009, Int. J. Bifurcation and Chaos, 

19, 2823 

Ozorio de Almeida, A. M. 1990, Hamiltonian Systems: Chaos and 

Quantization (Cambridge: Cambridge University Press) 
Patience, J., White, R. J., Ghez, A. M., et al. 2002, ApJ, 581, 654 
Pilat-Lohinger, E., & Dvorak, R. 2002, Celest. Mech. Dyn. Astron., 
82, 143 

Rabl, G., & Dvorak, R. 1988, A&A, 191, 385 

Raghavan, D., Henry, T. J., Mason, B. D., et al. 2006, ApJ, 646, 523 
Roy, A. E. 2005, Orbital Motion (Bristol and Philadelphia: Institute 
of Physics Publ.) 

Smith, R. H., & Szebehely, V. 1993, Celest. Mech. Dyn. Astron., 56, 
409 

Szebehely, V. 1967, Theory of Orbits (New York and London: 
Academic Press) 

Szenkovits, F., & Mako, Z. 2008, Celest. Mech. Dyn. Astron., 101, 273 
Tsonis, A. A. 1992, Chaos. From Theory to Applications (New York: 

Plenum Press) 
Yoshida, H. 1990 Phys. Lett. A, 150, 262 

Wolf, A., Swift, J. B., Swinney, H. L., & Vastano, J. A. 1985, Physica 
D, 16, 285 



Quarlcs ct al.: The instability transition for the restricted 3-body problem. III. 




D) 

5-1- 

v — ' 1 

a 

LU * 

1-6' 
a 

«_ 7 . 



^2 *3 



I I I I I I I I I I I I I I I I I I I I I 



8 10 



Time (binary orbits x 10 ) 




i 1 1 1 1 1 1 1 1 1 1 1 

12 3 4 
Time (binary orbits) 



Fig. 2. Case study for the initial planetary distance ratio po = 0.355 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p, = 0.3. 
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Fig. 3. Case study for the initial planetary distance ratio po = 0.474 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p — 0.3. 




■< 1 

o> 
o 

~ 



c 

CD 

c . 

o 

Q. 

X 

LU • 

> 
O 

c 

13 
Q. 
CO 
>, . 



b 




1 „ . ■ • i 1 i ■ i 1 . 







50 100 150 
Time (binary orbits) 




1 2 3 4 5 
Time (binary orbits) 



Fig. 4. Case study for the initial planetary distance ratio po = 0.595 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p — 0.3. 
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Fig. 5. Case study for the initial planetary distance ratio po = 0.290 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p = 0.5. 
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Fig. 6. Case study for the initial planetary distance ratio po = 0.370 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p — 0.5. 
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Fig. 7. Case study for the initial planetary distance ratio po = 0.400 with the planetary orbit in the synodic coordinate 
system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is p — 0.5. 
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Maximum Lyapunov Exponent t=10 5 yrs Maximum Lyapunov Exponent x=10 5 yrs 




Fig. 8. We depict the maximum Lyapunov exponent (colour coded) for various mass ratios [i and initial conditions po 
using a linear scale (left) or a logarithmic scale (right) where — n corresponds to 10 _n . The crosses denote cases where 
the simulation was terminated due to the planet being captured by one of the stars or being ejected from the system. The 
green, red, and black curves represent the initial conditions that cause the ZVC to open at LI, L2, and L3, respectively. 



